This file was last update on 2020-10-02
No. of visits to each transect/site.
abn%>%
group_by(Study, Taxa, Transect)%>%
summarise(n = length(unique(Date)))%>%
spread(Transect, n)
| Study | Taxa | T1A | T1B | T2 | T3 | T4 | T5 | T6 |
|---|---|---|---|---|---|---|---|---|
| 1998 - 2000 | Birds | 19 | 18 | 7 | 10 | 4 | 1 | 1 |
| 1998 - 2000 | Butterflies | 14 | 7 | 4 | 9 | 4 | 1 | 1 |
| 2016 - 2017 | Birds | 23 | 23 | 22 | 19 | 16 | 22 | 17 |
| 2016 - 2017 | Butterflies | 23 | 20 | 20 | 22 | 18 | 21 | 15 |
No. of individuals sampled
abn%>%
group_by(Study, Taxa, Transect)%>%
summarise(n = sum(Abundance))%>%
spread(Transect, n)
| Study | Taxa | T1A | T1B | T2 | T3 | T4 | T5 | T6 |
|---|---|---|---|---|---|---|---|---|
| 1998 - 2000 | Birds | 255 | 165 | 103 | 213 | 236 | 7 | 28 |
| 1998 - 2000 | Butterflies | 93 | 59 | 12 | 180 | 157 | 10 | 4 |
| 2016 - 2017 | Birds | 480 | 393 | 358 | 289 | 165 | 203 | 133 |
| 2016 - 2017 | Butterflies | 432 | 188 | 207 | 652 | 187 | 224 | 124 |
Uneveness in sampling effort across sites in studies may bias estimates of change.
library(vegan)
# input for vegan
veg_birds <- abn%>%
filter(Taxa == "Birds")%>%
group_by(Study, Transect, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = paste(Study, "_", Transect))%>%
group_by(sample)%>%
#mutate(N = sum(n, na.rm = T))%>%
#filter(N > 0)%>%#removing samples where no species were detected
dplyr::select(sample, Specific.Epithet, n)%>%
ungroup()%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
veg_butterfly <- abn%>%
filter(Taxa == "Butterflies")%>%
group_by(Study, Transect, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = paste(Study, "_", Transect))%>%
group_by(sample)%>%
#mutate(N = sum(n, na.rm = T))%>%
#filter(N > 0)%>%#removing samples where no species were detected
dplyr::select(sample, Specific.Epithet, n)%>%
ungroup()%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
## Species accumulation curves
specaccum(veg_butterfly)
## Species Accumulation Curve
## Accumulation method: exact
## Call: specaccum(comm = veg_butterfly)
##
##
## Sites 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000
## Richness 26.85714 39.51648 47.34615 52.96004 57.33017 60.91275 63.96037
## sd 12.83331 11.06085 9.33188 8.00650 7.02146 6.19857 5.49616
##
## Sites 8.00000 9.00000 10.00000 11.00000 12.00000 13.00000 14
## Richness 66.62438 68.99800 71.13986 73.08791 74.86813 76.50000 78
## sd 4.88039 4.31448 3.76625 3.21657 2.59888 1.76271 0
rareslope(veg_butterfly, sample = 200)
## 1998 - 2000 _ T1A 1998 - 2000 _ T1B 1998 - 2000 _ T2 1998 - 2000 _ T3
## 0.00000000 0.00000000 0.00000000 0.00000000
## 1998 - 2000 _ T4 1998 - 2000 _ T5 1998 - 2000 _ T6 2016 - 2017 _ T1A
## 0.00000000 0.00000000 0.00000000 0.04559341
## 2016 - 2017 _ T1B 2016 - 2017 _ T2 2016 - 2017 _ T3 2016 - 2017 _ T4
## 0.00000000 0.05079708 0.04937664 0.00000000
## 2016 - 2017 _ T5 2016 - 2017 _ T6
## 0.05173746 0.00000000
rarefy(veg_butterfly, sample = 200)
## 1998 - 2000 _ T1A 1998 - 2000 _ T1B 1998 - 2000 _ T2 1998 - 2000 _ T3
## 27.00000 22.00000 8.00000 37.00000
## 1998 - 2000 _ T4 1998 - 2000 _ T5 1998 - 2000 _ T6 2016 - 2017 _ T1A
## 24.00000 5.00000 3.00000 34.98930
## 2016 - 2017 _ T1B 2016 - 2017 _ T2 2016 - 2017 _ T3 2016 - 2017 _ T4
## 32.00000 35.65386 31.83323 31.00000
## 2016 - 2017 _ T5 2016 - 2017 _ T6
## 37.84749 27.00000
## attr(,"Subsample")
## [1] 200
rarebutts <- rarecurve(veg_butterfly)
abn%>%
group_by(Taxa, Study)%>%
summarise(richness = length(unique(Specific.Epithet)),
N = sum(Abundance, na.rm = T))
| Taxa | Study | richness | N |
|---|---|---|---|
| Birds | 1998 - 2000 | 70 | 1007 |
| Birds | 2016 - 2017 | 105 | 2021 |
| Butterflies | 1998 - 2000 | 45 | 515 |
| Butterflies | 2016 - 2017 | 66 | 2014 |
More species of both taxa were encountered in current study compared to the previous study. Difference in sampling effort is apparent accross two studies.
No. of species encountered across sites sampled in the two studies.
abn%>%
left_join(transects, by = c("Transect"))%>%
group_by(Taxa, Study, Transect)%>%
summarise(richness = length(unique(Specific.Epithet)))%>%
spread(Transect, richness)
| Taxa | Study | T1A | T1B | T2 | T3 | T4 | T5 | T6 |
|---|---|---|---|---|---|---|---|---|
| Birds | 1998 - 2000 | 47 | 37 | 28 | 36 | 29 | 4 | 11 |
| Birds | 2016 - 2017 | 66 | 59 | 53 | 47 | 31 | 38 | 28 |
| Butterflies | 1998 - 2000 | 27 | 22 | 8 | 37 | 24 | 5 | 3 |
| Butterflies | 2016 - 2017 | 41 | 32 | 36 | 44 | 31 | 39 | 27 |
Fewer species were encountered in sites with lower effort in the previous study, i.e., T5 and T6.
# Calculating diversity in each site over years
div_year <- abn%>%
group_by(Study, Taxa, Transect, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
spread(Specific.Epithet, n, fill = c(0))%>%
group_by(Study, Taxa, Transect)%>%
nest()%>%
mutate(H = map_dbl(data, ~vegetarian::d(.[-1], lev = "alpha", q = 1)))%>%
dplyr::select(-data)
# Summary
div_year%>%
group_by(Taxa, Study)%>%
skimr::skim(H)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Study, mean, sd)
Variable type: numeric
| Taxa | Study | mean | sd |
|---|---|---|---|
| Birds | 1998 - 2000 | 15.39 | 7.79 |
| Birds | 2016 - 2017 | 24.88 | 5.74 |
| Butterflies | 1998 - 2000 | 11.06 | 6.57 |
| Butterflies | 2016 - 2017 | 20.00 | 3.01 |
# t test
div_year%>%
group_by(Taxa)%>%
nest()%>%
mutate(test = map(data, ~t.test(H ~ Study, data = .)),
sumry = map(test, broom::tidy))%>%
dplyr::select(sumry)%>%
unnest()%>%
tstats()
| Taxa | 1998-2001 | 2016-2017 | parameter | statistic | p.value |
|---|---|---|---|---|---|
| Birds | 15.39287 | 24.87583 | 11.031229 | -2.593073 | 0.0249503 |
| Butterflies | 11.05543 | 19.99721 | 8.410531 | -3.273453 | 0.0105504 |
# plot
div_year%>%
ggplot(aes(Taxa, H, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5)+
labs(y = "Effective number of species")+
scale_fill_brewer(palette = "Greys")
ggsave(last_plot(), filename = "./Figures and Tables/figure2.png", height = 6, width = 8, dpi = 300)
The first order diversity of both taxa has increased significantly across the two surveys.
# sample info
samp_birds <- abn%>%
left_join(transects, by = c("Transect"))%>%
filter(Taxa == "Birds")%>%
dplyr::select(Study, Taxa, Transect, Habitat)%>%
distinct()%>%
mutate(sample = paste(Study, "_", Transect))
samp_butterfly <- abn%>%
left_join(transects, by = c("Transect"))%>%
filter(Taxa == "Butterflies")%>%
dplyr::select(Study, Taxa, Transect, Habitat)%>%
distinct()%>%
mutate(sample = paste(Study, "_", Transect))
# Similarity matrix by studies
sim_birds <- abn%>%
filter(Taxa == "Birds")%>%
group_by(Study, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = Study)%>%
ungroup()%>%
dplyr::select(sample, Specific.Epithet, n)%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
sim_butterfly <- abn%>%
filter(Taxa == "Butterflies")%>%
group_by(Study, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = Study)%>%
ungroup()%>%
dplyr::select(sample, Specific.Epithet, n)%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
# bray - curtis
data_frame(
Taxa = c("Birds", "Butterflies"),
Dissimilarity = c(as.numeric(vegdist(sim_birds, method = "bray")),
as.numeric(vegdist(sim_butterfly, method = "bray")))
)
| Taxa | Dissimilarity |
|---|---|
| Birds | 0.5495376 |
| Butterflies | 0.6710162 |
# ordination
NMDS_birds <- metaMDS(veg_birds, distance = "bray", k = 2, model = "local")
NMDS_butterfly <- metaMDS(veg_butterfly, distance = "bray", trymax = 1000, k = 2, model = "local")
# Extracting scores and adding habitat vars
ord_birds <- data.frame(scores(NMDS_birds))%>%
rownames_to_column("sample")%>%
left_join(samp_birds, "sample")
ord_butterfly <- data.frame(scores(NMDS_butterfly))%>%
rownames_to_column("sample")%>%
left_join(samp_butterfly, "sample")
ord <- bind_rows(ord_birds, ord_butterfly)
#Testing change in composition of birds
adonis(veg_birds ~ Study, data = samp_birds)
##
## Call:
## adonis(formula = veg_birds ~ Study, data = samp_birds)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Study 1 0.77547 0.77547 3.9632 0.24827 0.002 **
## Residuals 12 2.34798 0.19567 0.75173
## Total 13 3.12345 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Testing change in composition of butterflies
adonis(veg_butterfly ~ Study, data = samp_butterfly)
##
## Call:
## adonis(formula = veg_butterfly ~ Study, data = samp_butterfly)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Study 1 0.9191 0.91908 4.054 0.25253 0.001 ***
## Residuals 12 2.7205 0.22671 0.74747
## Total 13 3.6396 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plotting
ord%>%
ggplot(aes(NMDS1, NMDS2, col = Study, fill = Study))+
stat_ellipse(level = 0.95, geom = "polygon", linetype = "dashed", size = 1, alpha = 0.3)+
geom_point(aes(shape = Habitat), size = 3)+
guides(shape = guide_legend(nrow = 2, byrow = T))+
scale_color_brewer(palette = "Set1")+
scale_fill_brewer(palette = "Set1")+
facet_wrap(~Taxa)+
theme(legend.box = "vertical")
ggsave(last_plot(), filename = "./Figures and Tables/figure3.png", height = 6, width = 8, dpi = 300)
Both bird and butterfly communities compositions are significantly different compared to the previous study.
# compiling trophic guild data
host <- read.csv("./Data/Butterflies_host plant_280820.csv") #butterfly host plant data
diet <- read.csv("./Data/Birds_diet.csv") #bird diet data
trophic <- host%>%
full_join(diet)%>%
dplyr::select(Specific.Epithet, Guild)%>%
distinct()
# Calculating diversity of guilds
guilds <- abn%>%
left_join(trophic, by = c("Specific.Epithet"))%>%
filter(!Guild == "")%>%
group_by(Study, Taxa, Transect, Guild, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
spread(Specific.Epithet, n, fill = c(0))%>%
group_by(Study, Taxa, Transect, Guild)%>%
nest()%>%
mutate(H = map_dbl(data, ~vegetarian::d(.[-1], lev = "alpha", q = 1)))%>%
dplyr::select(-data)
# summary
guilds%>%
group_by(Taxa, Guild, Study)%>%
skimr::skim(H)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Guild, Study, mean, sd)
Variable type: numeric
| Taxa | Guild | Study | mean | sd |
|---|---|---|---|---|
| Birds | Carnivore | 1998 - 2000 | 3.01 | 1.14 |
| Birds | Carnivore | 2016 - 2017 | 4.67 | 2.62 |
| Birds | Frugivore | 1998 - 2000 | 3.37 | 1.51 |
| Birds | Frugivore | 2016 - 2017 | 4.21 | 0.60 |
| Birds | Grainivore | 1998 - 2000 | 2.44 | 1.43 |
| Birds | Grainivore | 2016 - 2017 | 2.83 | 1.09 |
| Birds | Insectivore | 1998 - 2000 | 7.74 | 3.06 |
| Birds | Insectivore | 2016 - 2017 | 13.58 | 4.82 |
| Birds | Omnivore | 1998 - 2000 | 3.90 | 1.78 |
| Birds | Omnivore | 2016 - 2017 | 5.88 | 1.15 |
| Butterflies | Generalist | 1998 - 2000 | 2.79 | 1.28 |
| Butterflies | Generalist | 2016 - 2017 | 5.62 | 1.91 |
| Butterflies | Grass Specialist | 1998 - 2000 | 1.87 | 0.57 |
| Butterflies | Grass Specialist | 2016 - 2017 | 2.98 | 0.74 |
| Butterflies | Herb Specialist | 1998 - 2000 | 2.73 | 0.93 |
| Butterflies | Herb Specialist | 2016 - 2017 | 4.09 | 0.86 |
| Butterflies | Liana Specialist | 1998 - 2000 | 1.00 | 0.00 |
| Butterflies | Shrub Specialist | 1998 - 2000 | 3.54 | 1.71 |
| Butterflies | Shrub Specialist | 2016 - 2017 | 4.78 | 1.50 |
| Butterflies | Tree Specialist | 1998 - 2000 | 4.63 | 3.02 |
| Butterflies | Tree Specialist | 2016 - 2017 | 6.88 | 1.63 |
# t - test
guilds%>%
filter(Guild != "Liana Specialist")%>% #removing due to lack of obs
complete(nesting(Study, Taxa), Transect, Guild, fill = list(H = 0))%>%
group_by(Taxa, Guild)%>%
nest()%>%
mutate(test = map(data, ~t.test(H ~ Study, data = .)),
sumry = map(test, broom::tidy))%>%
dplyr::select(sumry)%>%
unnest()%>%
tstats()
| Taxa | Guild | 1998-2001 | 2016-2017 | parameter | statistic | p.value |
|---|---|---|---|---|---|---|
| Birds | Carnivore | 3.013676 | 3.334886 | 40.52328 | -0.6191168 | 0.5393047 |
| Birds | Frugivore | 3.367593 | 4.211182 | 62.64226 | -3.8933034 | 0.0002426 |
| Birds | Grainivore | 2.440818 | 2.833390 | 89.94433 | -1.6340146 | 0.1057514 |
| Birds | Insectivore | 6.631103 | 13.578874 | 88.97385 | -7.9838426 | 0.0000000 |
| Birds | Omnivore | 3.901475 | 5.876942 | 82.45283 | -6.9810904 | 0.0000000 |
| Butterflies | Generalist | 1.990058 | 5.622563 | 77.56402 | -9.7435153 | 0.0000000 |
| Butterflies | Grass Specialist | 1.335180 | 2.975740 | 57.94419 | -8.6211871 | 0.0000000 |
| Butterflies | Herb Specialist | 1.560185 | 4.092249 | 35.98571 | -8.2392106 | 0.0000000 |
| Butterflies | Shrub Specialist | 3.542509 | 4.780687 | 94.34174 | -4.0786608 | 0.0000946 |
| Butterflies | Tree Specialist | 4.631978 | 6.879746 | 73.68039 | -4.8985262 | 0.0000056 |
# plot
guilds%>%
ggplot(aes(reorder(Guild, H), H, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5, position = position_dodge(preserve = "single"))+
labs(y = "Effective number of species", x = "Trophic Guild")+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
scale_fill_brewer(palette = "Greys")+
facet_wrap(~Taxa, scales = "free_x", ncol = 1)
ggsave(last_plot(), filename = "./Figures and Tables/figure4.png", height = 9, width = 8, dpi = 300)
Carnivorous, Grainivorous birds have reduced in the comunity, compared to insectivorous birds which have increased. Similarly, Shrub specialist and tree specialist butterflies have reduced where grass specialists and herb specialist have increased. Suprsingly, generalist species of both taxa appear to be unchanged.
-Can we incorporate abunndance into community specialisation index?
# trophic specialisation based on Juliard et al. 2006
butterflies_TSI <- host%>%
mutate(H = length(unique(Hostplant.Family)))%>%
group_by(Specific.Epithet)%>%
mutate(h = length(unique(Hostplant.Family)))%>%
ungroup()%>%
mutate(TSI = sqrt((H/h)-1))
birds_TSI <- diet%>%
mutate(H = length(unique(Prey)))%>%
group_by(Specific.Epithet)%>%
mutate(h = length(unique(Prey)))%>%
ungroup()%>%
mutate(TSI = sqrt((H/h)-1))
CSI.T <- abn%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
group_by(Study, Taxa, Transect)%>%
distinct()%>%
left_join(butterflies_TSI, by = c("Specific.Epithet"))%>%
left_join(birds_TSI, by = c("Specific.Epithet"))%>%
mutate(TSI = ifelse(is.na(TSI.x), TSI.y, TSI.x))%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet, TSI)%>%
summarise(CSI.T = mean(TSI, na.rm = T))
# summary
CSI.T%>%
group_by(Taxa, Study)%>%
skimr::skim(CSI.T)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Study, mean, sd)
Variable type: numeric
| Taxa | Study | mean | sd |
|---|---|---|---|
| Birds | 1998 - 2000 | 2.59 | 0.16 |
| Birds | 2016 - 2017 | 2.53 | 0.09 |
| Butterflies | 1998 - 2000 | 5.34 | 0.39 |
| Butterflies | 2016 - 2017 | 5.07 | 0.24 |
# t - test
CSI.T%>%
group_by(Taxa)%>%
nest()%>%
mutate(test = map(data, ~t.test(CSI.T ~ Study, data = .)),
sumry = map(test, broom::tidy))%>%
dplyr::select(sumry)%>%
unnest()%>%
tstats()
| Taxa | 1998-2001 | 2016-2017 | parameter | statistic | p.value |
|---|---|---|---|---|---|
| Birds | 2.588009 | 2.531490 | 9.462567 | 0.8103492 | 0.4376441 |
| Butterflies | 5.338563 | 5.065191 | 10.056310 | 1.5727722 | 0.1466765 |
CSI.T%>%
ggplot(aes(Taxa, CSI.T, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5)+
labs(y = "Community specialisation index \n (Trophic)")+
scale_fill_brewer(palette = "Greys")+
scale_y_sqrt()
Niche width of the community did not change over the two studies.
HSI <- abn%>%
left_join(transects, by = c("Transect"))%>%
group_by(Study, Taxa, Habitat, Specific.Epithet)%>%
# calculating abundance of each species in each habitat type
summarise(n = sum(Abundance, na.rm = T))%>%
group_by(Study, Taxa)%>%
mutate(N = sum(n), # total number of indviduals encountered
rel.prop = n/N)%>% # relative proportion of species
group_by(Study, Specific.Epithet)%>%
summarise(HSI = sd(rel.prop)/mean(rel.prop))%>%
replace_na(replace = list(HSI = 0))
CSI.H <- abn%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
group_by(Study, Taxa, Transect)%>%
distinct()%>%
left_join(HSI, by = c("Study", "Specific.Epithet"))%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet, HSI)%>%
summarise(CSI.H = mean(HSI, na.rm = T))
# summary
CSI.H%>%
group_by(Taxa, Study)%>%
skimr::skim(CSI.H)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Study, mean, sd)
Variable type: numeric
| Taxa | Study | mean | sd |
|---|---|---|---|
| Birds | 1998 - 2000 | 0.54 | 0.10 |
| Birds | 2016 - 2017 | 0.49 | 0.02 |
| Butterflies | 1998 - 2000 | 0.61 | 0.09 |
| Butterflies | 2016 - 2017 | 0.57 | 0.03 |
# t - test
CSI.H%>%
group_by(Taxa)%>%
nest()%>%
mutate(test = map(data, ~t.test(CSI.H ~ Study, data = .)),
sumry = map(test, broom::tidy))%>%
dplyr::select(sumry)%>%
unnest()%>%
tstats()
| Taxa | 1998-2001 | 2016-2017 | parameter | statistic | p.value |
|---|---|---|---|---|---|
| Birds | 0.5363977 | 0.4942803 | 6.719090 | 1.068760 | 0.3220632 |
| Butterflies | 0.6110527 | 0.5692785 | 7.782004 | 1.210166 | 0.2616879 |
CSI.H%>%
ggplot(aes(Taxa, CSI.H, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5)+
labs(y = "Community specialisation index \n (Habitat)")+
scale_fill_brewer(palette = "Greys")
Habitat niche of the community did not change over the study periods.
season <- abn%>%
mutate(Month = month(dmy(Date), label = T))%>%
group_by(Study, Month, Taxa)%>%
summarise(n = sum(Abundance, na.rm=T))%>%
group_by(Study, Taxa)%>%
mutate(N = sum(n),
rel.prop = n/N)
season%>%
ggplot(aes(Month, rel.prop, col = Taxa, shape = Taxa))+
geom_point(size = 3)+
geom_path(aes(group = Taxa), size = 1, linetype = "dashed")+
scale_color_brewer(palette = "Set1")+
facet_wrap(~Study, ncol = 1)
library(raster)
library(sf)
# Raw point data
points <- read.csv("./Data/sites.csv")%>%
drop_na()%>%
dplyr::select(-Location)
# Converting to SF points
points_sf <- st_as_sf(points, coords = c('Lon', 'Lat'))
st_crs(points_sf) <- 4326
# Making transect lines
lines <- points_sf%>%
mutate(Transect = as.factor(Transect))%>%
group_by(Transect)%>%
summarise()%>%
st_cast("LINESTRING")
# getting study extent
study_ext <- st_bbox(points_sf)
ext <- st_as_sfc(study_ext)
# getting osm map
library(osmdata)
# Tamhini WLS
tamhini_raw <- opq(bbox = study_ext, timeout = 100)%>%
add_osm_feature(key = 'boundary', value = 'protected_area')%>%
osmdata_sf()%>%
unique_osmdata()
tamhini <- tamhini_raw$osm_multipolygons
st_crs(tamhini) <- 4326
# Roads
tamhini_area <- st_bbox(tamhini)
roads_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
add_osm_feature(key = 'highway')%>%
osmdata_sf()%>%
unique_osmdata()
roads <- roads_raw$osm_lines
st_crs(roads) <- 4326
# Landuse
landuse_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
add_osm_feature(key = 'natural')%>%
osmdata_sf()%>%
unique_osmdata()
forest <- landuse_raw$osm_polygons%>%
filter(natural == "wood")
water <- landuse_raw$osm_multipolygons
water_small <- landuse_raw$osm_polygons%>%
filter(natural == "water")
st_crs(forest) <- 4326
# waterways
waterway_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
add_osm_feature(key = 'waterway')%>%
osmdata_sf()%>%
unique_osmdata()
waterways <- waterway_raw$osm_lines
st_crs(waterways) <- 4326
# elevation
library(elevatr)
tam_elev <- elevatr::get_elev_raster(locations = tamhini, z = 8)
tam_contours <- raster::rasterToContour(tam_elev, levels = seq(0, 1000, 100))
# making final map
library(gridExtra)
fig1 <- grid.arrange(fig1a, fig1b, fig1c, legend, layout_matrix = rbind(c(1, 1, 1),
c(1, 1, 1),
c(2, 3, 4)))
ggsave(fig1, filename = "./Figures and Tables/figure1.png", height = 8, width = 8)